df <- read.csv2(file = "data/ZNOandVoating/input_2016.csv", header = 1, sep = ",", dec = ".", fileEncoding = "CP1251", stringsAsFactors = FALSE)
names(df) <- c("regname", "ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")

head(df, 3)

##             regname ukr math  pro_ukr  radical opposition   small
## 1 Вінницька область 167  128 43.00824 6.005472   1.312608 7.75368
## 2 Вінницька область 137    0 43.00824 6.005472   1.312608 7.75368
## 3 Вінницька область 141    0 43.00824 6.005472   1.312608 7.75368
##   not_voted
## 1     41.92
## 2     41.92
## 3     41.92

n_for_delete <- nrow(df[df$regname == "Донецька область" | df$regname == "Луганська область", ])
print(as.character(round(n_for_delete * 100 / nrow(df), 2) %&% "% записів, що відповідають Донецькій та Луганській облатям"))

## [1] "5.87% записів, що відповідають Донецькій та Луганській облатям"

df <- df[df$regname != "Донецька область" & df$regname != "Луганська область", ]

ЗНО з математики та української мови

Виведемо загальні характеристики даних ЗНО з математики та української мови:

describe(df[, c("math", "ukr")])

## df[, c("math", "ukr")] 
## 
##  2  Variables      252262  Observations
## ---------------------------------------------------------------------------
## math 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       55     0.77    53.53    69.96        0        0 
##      .25      .50      .75      .90      .95 
##        0        0      123      157      172 
## 
## lowest :   0 100 104 107 111, highest: 196 197 198 199 200
## ---------------------------------------------------------------------------
## ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       83    0.997    125.4    58.72        0        0 
##      .25      .50      .75      .90      .95 
##      111      137      165      184      189 
## 
## lowest :   0.0 100.0 101.0 103.0 104.0, highest: 197.0 198.0 199.0 199.5 200.0
## ---------------------------------------------------------------------------

Із таблиці вище видно, що в записах нема пропущених значень. Судячи із значень квантилів, половина значень ЗНО з математики є нулями. Варто відмітити, що значно менша частка значень 0 у ЗНО з української мови (близько 10%).

Подивимося на розподіли:

p1 <- plot_ly(x = df[, "math"], type = "histogram", name = "math")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram", name = "ukr")

plotly::subplot(p1, p2) %>% layout(title = "EIT: Ukranian language and mathematics")

remove(p1); remove(p2)

Розподіли результатів ЗНО з української та математики значно різняться.

Подивимося на співвідношення людей, котрі отримали 0 по математиці та українській одночасно

print("Не подолали поріг або не здавали ЗНО відносно всіх абітурієнтів: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df), 2)) %&% "%")

## [1] "Не подолали поріг або не здавали ЗНО відносно всіх абітурієнтів: 12.72%"

print("Не подолали поріг або не здавали ЗНО відносно тих, хто не склав математику: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df[(df$math == 0), ]), 2)) %&% "%")

## [1] "Не подолали поріг або не здавали ЗНО відносно тих, хто не склав математику: 20.77%"

print("Не подолали поріг або не здавали ЗНО відносно тих, хто не склав українську: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df[(df$ukr == 0), ]), 2)) %&% "%")

## [1] "Не подолали поріг або не здавали ЗНО відносно тих, хто не склав українську: 91.19%"

Із останнього напрошуєтся висновок, що ті, хто не склали іспит із української, скоріш за все не склали й математику. Що є надзвичайно логічно: для вступу в університет абітурієнту необхідно отримати певний мінімум балів із кожного ЗНО. Першою завжди здається українська. Тому, не склавши українську, сенс складати математику залишається у винятковим.

Поглянемо на кореляції між результатами ЗНО, політичними силами:

corrolation_m <- cor(df[, -1], method = "pearson")
plot_ly(x = rownames(corrolation_m), y = colnames(corrolation_m), z = corrolation_m, colors = "Greys", type = "heatmap") %>%
  layout(title = "Pearson Correlation")

# remove(corrolation_m)

Відмітимо майже 0 кореляцію між результатами ЗНО з української мови та про-українським політичним напрямком та слабку кореляцію між результатми математики та української.

Досить висока кореляція між результатми не голосування та голосуванням за малі партії. Така ситуація може пояснюватися тим, що населення не підтримує основних лідерів, а тому не йде голосувати або голосує спеціально за тих, хто не пройде.

Видалимо стрічки, що відповідають спостереженням із 0 хочаб в одному ЗНО та подивимося на розподіл.

df <- df[(df$math != 0) & (df$ukr != 0) ,]
p1 <- plot_ly(x = df[, "math"], type = "histogram", name = "math")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram", name = "ukr")

plotly::subplot(p1, p2) %>% layout(title = "Distribution of the passed EIT",
  annotations = list(list(x = 0.2, y = 1.05, text = "math", showarrow = F, xref = 'paper', yref = 'paper'),
                     list(x = 0.8, y = 1.05, text = "ukr", showarrow = F, xref = 'paper', yref = 'paper')))

remove(p1); remove(p2);

Відмітимо, що наразі дуже відсіялася частка людей, котрі погано здали українську мову (було в околі 3-х тисяч, стало в около 700).

Гістограма розсіювання, в свою чергу, виглядає наступним чином:

plot_ly(data = df, x = ~math, y = ~ukr, mode = "markers", type = "scatter", marker=list(opacity = 0.7)) %>% layout(title = "Passed EIT results")

Варто відмітити, що для абітурієнтів, що склали математику краще ніж 160 балів, помітна тенденція, як правило, кращого складання ЗНО з математики.

Подивимося на розподіли по кожній політичній силі окремо (якщо конкретно обрана політична сила має для даного спостереження найбільший вплив (найбільшу ймовірність), то дане спостереження відносимо до ціє політичної сили)

• математика

W <- as.data.frame(df[, -c(1, 2, 3)])
W[, "max_value"] <- apply(W, 1, max)

all_plots <- list()
annotations_list <- list()
names_ <- names(W[, -ncol(W)])

for (i in seq_along(names_)) {
  name <- names_[i]
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {
    all_plots[[name]] <- plot_ly(x = d[, "math"], type = "histogram", name = name)
    annotations_list[[i]] <- list(x = 0.15*i, y = 1.005, text = name, showarrow = F, xref = 'paper', yref = 'paper')
  }
}

plotly::subplot(all_plots) %>%
  layout(title = "Mathematics EIT according to political strength", annotations = annotations_list)

remove(d); remove(all_plots); remove(annotations_list); remove(i); remove(name)

Відмітимо, що візуально розподіли балів виглядають майже однаково для про-украхнських регіонів та тих, де в основному, не голосували.

• українська мова

all_plots <- list()
annotations_list <- list()

for (i in seq_along(names_)) {
  name <- names_[i]
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {
    all_plots[[name]] <- plot_ly(x = d[, "ukr"], type = "histogram", name = name)
    annotations_list[[i]] <- list(x = 0.15*i, y = 1.005, text = name, showarrow = F, xref = 'paper', yref = 'paper')
  }

}

plotly::subplot(all_plots[["pro_ukr"]], all_plots[["not_voted"]]) %>%
  layout(title = "Ukr EIT according to political strength", annotations = annotations_list)

remove(d); remove(all_plots); remove(annotations_list); remove(names_); remove(name); remove(i)

Розподіли результатів ЗНО з уркаїнської візуально більш різняця для різних політичних вподобань.

for (name in names(W[, -ncol(W)])) {
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {WVPlots::ScatterHist(d, "math", "ukr", name, smoothmethod = "none", estimate_sig = FALSE, minimal_labels = TRUE)}
}

remove(d); remove(W); remove(name); remove(df)

Із останньої пари зображень здається, що нижній поріг результатів української та математики (для абітурієнтів із математикою > 160 балів) вищий у випадку про-українських регіонів. Така ситуація можлива також через різну кількість абітурієнтів у порівнювальних регіонах.

Результати голосування

Виведемо загальні характеристики кожної політичної сили:

df <- read.csv2(file = "data/ZNOandVoating/input_2016.csv", header = 1, sep = ",", dec = ".", fileEncoding = "CP1251", stringsAsFactors = FALSE)
names(df) <- c("regname", "ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")
df <- df[df$regname != "Донецька область" & df$regname != "Луганська область", ]
describe(df[, -(1:3)])

## df[, -(1:3)] 
## 
##  5  Variables      252262  Observations
## ---------------------------------------------------------------------------
## pro_ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       23    0.997    33.51    13.31    15.71    16.23 
##      .25      .50      .75      .90      .95 
##    21.52    34.08    42.46    48.82    53.87 
## 
## lowest : 15.71244 16.22691 19.37149 20.13312 21.00261
## highest: 43.00824 44.77893 48.82355 50.41795 53.87200
## ---------------------------------------------------------------------------
## radical 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       23    0.997    6.703    2.561    3.047    3.832 
##      .25      .50      .75      .90      .95 
##    4.992    6.580    8.477    8.912   10.428 
## 
## lowest :  3.046992  3.831828  3.834072  4.523960  4.619912
## highest:  8.886996  8.911950 10.084956 10.427880 11.555155
## ---------------------------------------------------------------------------
## opposition 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       23    0.997    4.295    4.541   0.4165   0.4970 
##      .25      .50      .75      .90      .95 
##   1.3126   2.0612   7.1334  11.6156  14.5749 
## 
## lowest :  0.344142  0.416508  0.497000  0.953295  1.042825
## highest:  6.796640  7.133360 10.466742 11.615622 14.574912
## ---------------------------------------------------------------------------
## small 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       23    0.997    9.319    2.148    6.779    7.189 
##      .25      .50      .75      .90      .95 
##    7.754    9.167    9.966   13.113   13.113 
## 
## lowest :  5.856787  6.778902  7.189344  7.360584  7.496422
## highest: 10.890088 11.198572 11.346280 13.112736 13.519935
## ---------------------------------------------------------------------------
## not_voted 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   252262        0       23    0.997    46.17    9.604    30.00    35.15 
##      .25      .50      .75      .90      .95 
##    41.92    44.14    52.81    57.20    60.48 
## 
## lowest : 30.00 31.72 35.15 36.27 39.79, highest: 54.68 55.32 57.20 58.64 60.48
## ---------------------------------------------------------------------------

Одразу видно, що 23 унікальні значення в кожному стовпчику відповідають 22 областям України та м. Києву. Варто відмітити, що відсутні області, із переважними радикальними або опозиційними силами. Також немає регіонів, у яких голоса віддані виключно за малі політичні партії.

Розподіл абітурієнтів по областях:

ab <- plyr::ddply(.data = df, .variables = .(pro_ukr, radical, opposition, small, not_voted), .fun = nrow)
plot_ly(ab, y = ~V1, type = "bar") %>%
  layout(title = "Abiturients distribution according regions",
         xaxis = list(title = "Region ID"), yaxis = list(title = ""))

Варто відмітити, що кілька регіонів надто відрізняються чисельністю. Через таку зміщенність даних можуть бути зміщені результати.

Із розподілу відсотків підтримки кожної політичної сили відповідно до регіону видно, що більшість регіонів є політично неактивними. Також варто відмітити, що розподіли голосів малих партії, у більшості випадків, збігається з розподілом голосів радикальної партії.

plot_ly(y = ab[, "pro_ukr"], type = "bar", name = "pro_ukr") %>%
  add_trace(y = ab[, "not_voted"], name = "not_voted") %>%
  add_trace(y = ab[, "radical"], name = "radical") %>%
  add_trace(y = ab[, "opposition"], name = "opposition") %>%
  add_trace(y = ab[, "small"], name = "small")

remove(ab); remove(df)

Висновок аналізу даних

На основі порівняння розподілів можна зробити висновок про те, що дані можуть бути описані за допомогою моделі суміші зі змінними концентраціями. Проте, зважаючи на кореляції та відсоткові співвідношення, можливо, модель суміші зі змінними концентраціями, буде не надто ефективною.

Визначення параметрів згладжування

Скористаємося 5-ти фолдовою кросс-валідацією для знаходження оптимальних параметрів згладжування для кожної моделі (проведемо описаний нижче есперимент 5 разів).

Для визначення оптимального параметра використаємо модуль МНК зважений ваговими коефіцієнтами:

\[MSE^{(m)} = \frac{1}{N} | \sum_{i=1}^{N}{a_{i:N}^{(m)}(Y_i - \hat{g}^{(m)}(X_i))^2} |\] Зауважимо, що модуль використовується для того щоб підчас підрахунків середнього зваженого MSE по всіх фолдах фіксованої компоненти, фіксованого h, не отримати 0 (який може бути хибно протрактований).

Отже, визначення оптимальних параметрів h відбувається за схемою описаною нижче.

• Зафіксуємо h, порахуємо зважений МНК для кожного фолда та знайдемо середнє та дисперсію по 5-ти фолдам для кожної компоненти. (Так прорахуємо для кожного фіксованого h.)

• Відкинемо 25% значень середніх, які відповідають найбльшим значенням дисперсії. Таким чином залишаємо лише ті результати, на яких моделі показали більш стабільну поведінку.

• Серед значень середніх результатів моделей визначаємо h, які відповідають найменшим середнім зваженим МНК (окремо для кожної компоненти).

# h_range <- c(seq(0.1, 1, 0.2), seq(1, 5, 0.5), 8, 10)
# 
# for (i in 1:5) {
#   print("***** " %&% as.character(i) %&% " ******")
#   set.seed(i)
#   cv_split <- cross_validation_split(train, k = 5)
#   GNW_cv_results <- GNWcv_across_h(h_range = h_range,
#                                  cv_df_split = cv_split,
#                                  X_colname = "math", Y_colname = "ukr",
#                                  W_colname = c("pro_ukr", "radical", "opposition", "small", "not_voted"),
#                                  use_parallel = TRUE)
#   write.csv(GNW_cv_results, "data/computation_results/2016cc_data/5cv_mean_std_" %&% as.character(i) %&% ".csv")
# }
# 
# # h_opt <- optimal_h(GNW_cv_results)
# remove(h_range); remove(i); remove(cv_split); remove(GNW_cv_results)

Результат середніх та диспресії збережемно в змінній \(GNW\_cv\_results\), а оптимальні \(h\), серед запропонованих для перебору, у змінну \(h\_opt\).

data_path <- "data/computation_results/2016cc_data/"
file_names <- dir(data_path) #where you have your files
file_names <- data_path %&% file_names

cv_results <- lapply(file_names, 
                     function(x){
                       df <- as.matrix(read.csv(x, row.names = 1))
                       l <- list(df=df , h=optimal_h(df))
                       return(l)
                       })

cv_results_h <- lapply(cv_results, function(x){x$h})
cv_results_mean_std <- lapply(cv_results, function(x){x$df})

print(cv_results_h)

## [[1]]
##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 438.4071 454.3044 1636.697 1035.322 351.3797
## [2,]   2.5000  10.0000    0.100    5.000   2.5000
## 
## [[2]]
##          [,1]     [,2]    [,3]     [,4]     [,5]
## [1,] 752.0023 610.5381 1900.39 522.0256 1014.304
## [2,]   5.0000   5.0000   10.00   8.0000    4.000
## 
## [[3]]
##         [,1]     [,2]    [,3]    [,4]     [,5]
## [1,] 636.954 574.1809 770.329 409.463 539.1719
## [2,]   4.500  10.0000   0.100   3.500   0.1000
## 
## [[4]]
##         [,1]     [,2]    [,3]     [,4]     [,5]
## [1,] 261.346 488.8267 1669.38 476.3409 145.2318
## [2,]   3.000   5.0000    2.50   4.5000   4.0000
## 
## [[5]]
##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 485.6412 512.2064 1866.826 271.5143 205.7276
## [2,]   3.5000   5.0000    5.000   8.0000   2.5000

Порівняння результатів кросс-валідації та прогнозування на відкладеній вибірці

Порівняємо тепер результати кросс-валідації з результатми прогнозування на тестовій частині даних.

# train <- read.csv("data/computation_results/train.csv", row.names = 1, stringsAsFactors = FALSE)
# test <- read.csv("data/computation_results/test.csv", row.names = 1, stringsAsFactors = FALSE)
# h_opt <- cv_results_h[[4]][2, ]
# 
# GNW <- GeneralisedNadarayaWatson$new()
# GNW$train(X_train = as.numeric(train[, "math"]),
#           Y_train = as.numeric(train[, "ukr"]),
#           W_train = as.matrix(train[, -c(1, 2)]))
# GNW$run_cluster()
# prediction_with_acoeff <- GNW$predict_in_parallel(X_test = as.numeric(test[, "math"]),
#                                                   W_test = as.matrix(test[, -c(1, 2)]), h = h_opt)
# GNW$stop_cluster()
# 
# prediction <- prediction_with_acoeff[["prediction"]]
# A_coeff <- prediction_with_acoeff[["A_test"]]
# 
# write.csv(prediction, "data/computation_results/prediction2016v.csv")
# write.csv(A_coeff, "data/computation_results/A_coeff2016v.csv")
# 
# remove(GNW); remove(train)

Порахуємо зважену суму МНК для прогнозу зробленого на тестових даних.

train <- read.csv("data/computation_results/train.csv", row.names = 1, stringsAsFactors = FALSE)
test <- read.csv("data/computation_results/test.csv", row.names = 1, stringsAsFactors = FALSE)

prediction <- as.matrix(read.csv("data/computation_results/prediction2016v.csv", row.names = 1, stringsAsFactors = FALSE))
A_coeff <- as.matrix(read.csv("data/computation_results/A_coeff2016v.csv", row.names = 1, stringsAsFactors = FALSE))

WMSE <- weighted_MSE(Y_true = as.numeric(test[, "ukr"]), Y_predicted = prediction, A_coeff = A_coeff)
names(WMSE) <- rep("comp.") %&% as.character(1:5)
WMSE

##      comp.1      comp.2      comp.3      comp.4      comp.5 
##    517.8019   -449.8565 366616.4095   -553.3369  -1092.8865

Відмітимо, що найменше значення зваженого МНК для моделей, що відповідають проукраїнській силі. Нагадаємо, найбільша по модулю кореляція спостерігалася між результатами ЗНО з української мовои та проукр. силами й тими, хто не проголосував. Також відмітимо, що раніше ми змогли виділити лише регіони, де переважає проукр напрямок та ті, в яких переважає не голосування. Тобто не виявилось регіонів, у яких переважають голоси за радикалів, за опозицію або за малі партії.

Порівняємо зважені МНК результатів прогнозування з результатами крос-валідації, використовуючи відповідними параметрами згладжування.

h_opt <- cv_results_h[[4]][2, ]
GNW_cv_results <-cv_results_mean_std[[4]]

cv_mean_results <- sapply(1:length(h_opt),
                          function(i){
                            GNW_cv_results[which("mean_" %&% as.character(h_opt[i]) == rownames(GNW_cv_results)), i]
                            })
cv_test_comparing <- rbind(cv_mean_results, WMSE)
rownames(cv_test_comparing) <- c("cross-validation", "test")
colnames(cv_test_comparing) <- 1:5
print(cv_test_comparing)

##                         1         2         3         4          5
## cross-validation 261.3460  488.8267   1669.38  476.3409   145.2318
## test             517.8019 -449.8565 366616.41 -553.3369 -1092.8865

Із останньої таблиці видно, що квадрати залишків першої моделі, порахованих для кросс-валідації та на відкладеній тестовій вибірці не надто різняця у порівнянні з іншими компонентами.

Це свідчить про більш-менш вдалий підбір параметрів згладжування для моделі, що відповідає проукраїнському настрою.

Для 3 інших моделей результати - невтішні. Останнє може бути зумовлене слабкою представленістю відповідних класів у вибірці. Можливо є сенс розглянути 3-х компонентну суміш: проукраїнську направленість, “не проукраїнську” (об’єднати результати голосуваня за малі партії, опозицію, радикалів) та “байдужих”.

Візуалізація результатів

Зобразмо спочатку загальну картину: гістограму розсіювання з усіма прознозами для всіх моделей:

plot_ly(data = test, x = ~math, y = ~ukr, mode = "markers", type = "scatter", name="EIT results") %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 1], name = "pro_ukr",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 2], name = "radical",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 3], name = "opposition",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 4], name = "small",  line=list(color="yellow"), marker=list(color="yellow", line=list(color="yellow")), mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 5], name = "not_voted",  mode = 'lines+markers') %>%
  layout(title = "Passed EIT results (test dataset)")

Окремо зобразимо для абітурієнтів, які проживають в областях, що ймовірніше наслідують проукраїнську групу:

train <- read.csv("data/computation_results/train.csv", row.names = 1, stringsAsFactors = FALSE)
test <- read.csv("data/computation_results/test.csv", row.names = 1, stringsAsFactors = FALSE)

prediction <- as.matrix(read.csv("data/computation_results/prediction2016v.csv", row.names = 1, stringsAsFactors = FALSE))
A_coeff <- as.matrix(read.csv("data/computation_results/A_coeff2016v.csv", row.names = 1, stringsAsFactors = FALSE))

Wtest <- as.data.frame(cbind(prediction, test[, -c(1, 2)]))
Wtest[, "max_value"] <- apply(Wtest[, -(1:5)], 1, max)

pro_ukr <- test[Wtest[, "pro_ukr"] == Wtest["max_value"], ]
p1 <- plot_ly(data = pro_ukr, x = ~math, y = ~ukr, mode = "markers", type = "scatter", name = "EIT results in pro_ukr regions") %>%
  add_trace(x = test[rownames(pro_ukr), "math"],
            y = Wtest[rownames(pro_ukr), 1], name = "pro_ukr model",  mode = 'markers')
p1

Відмітимо, що основні тренди моделі для проукраїнських регіонів не надто збереглися.